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ABSTRACT 

The propagation of cosmological ionization fronts during the reionization of the universe is strongly influenced 
by small-scale gas inhomogeneities due to structure formation. These inhomogeneities include both collapsed 
minihalos, which are generally self-shielding, and lower-density structures, which are not. The minihalos are 
dense and sufficiently optically-thick to trap intergalactic ionization fronts, blocking their path and robbing them 
of ionizing photons until the minihalo gas is expelled as an evaporative wind. The lower-density structures do not 
trap these fronts, but they can slow them down by increasing the overall recombination rate in the intergalactic 
medium (IGM). In this paper we study the effects of both types of inhomogeneities, including nonlinear clustering 
effects, and we find that both IGM clumping and collapsed minihalos have significant yet qualitatively different 
impacts on reionization. While the number density of minihalos on average increases strongly with time, the 
density of minihalos inside H II regions around ionizing sources is largely constant. Thus the impact of minihalos 
is essentially to decrease the number of ionizing photons available to the IGM at all epochs, which is equivalent 
to a reduction in the luminosity of each source. On the other hand, the effect of IGM clumping increases strongly 
with time, slowing down reionization and extending it. Thus while the impact of minihalos is largely degenerate 
with the unknown source efficiency, IGM clumping can help significantly in reconciling the recent observations 
of cosmic microwave background polarization with quasar absorption spectra at z ~ 6, which together point to an 
early but extended reionization epoch. 

Subject headings: hydrodynamics — radiative transfer — galaxies: halos — galaxies: high-redshift — intergalactic 
medium — cosmology: theory 



1. INTRODUCTION 

Recent polarization observations of the cosmic microwave 
background by the Wilkinson Microwave Anisotropy Probe 
(WMAP) imply that reionization was fairly advanced at z K ~ 15 
(Kogut et al. 2003). This came as a surprise. The prior detec- 
tion of the Gunn-Peterson effect in the spectra of high-redshift 
quasars had suggested that reionization was only just ending at 
z — 6 (White et al 2003; Fan et al. 2004). That was consistent 
with predictions of the most accurate numerical simulations in 
the current ACDM paradigm, which had all predicted this tran- 
sition at z K £8-10 (Ciardi et al. 2000; Gnedin 2000a; Ra- 
zoumov et al. 2002; Ciardi et al. 2003). Despite many poorly 
understood details concerning the star formation rate, the es- 
cape fraction of ionizing radiation, and the differences in nu- 
merical treatments of reionization, z re ~ 15 had seemed un- 
likely, and such an extended period of reionization, impossible. 

Now the race is on to reconcile the early onset of reionization 
suggested by WMAP with the high-redshift Gunn-Peterson ef- 
fect, which implies neighboring ionized patches finally grew to 
overlap at z ~ 6 (Haiman & Holder 2003; Cen 2003; Wyithe & 
Loeb 2003; Ciardi et al. 2003). One suggestion is that the uni- 
verse had two reionization epochs but recombined in between 
(Cen 2003), yet this ignores the unavoidable spread in redshifts 
intrinsic to any such IGM transition (Scannapieco, Schneider, 
& Ferrara 2003; Barkana & Loeb 2004; Furlanetto & Loeb 
2005). Other suggestions involve fine-tuning the ionizing pho- 
ton emissivity for different source halo masses, the escape frac- 
tion, and the (possibly metalicity-dependent) Initial Mass Func- 
tion (IMF), in ways intended to accelerate early ionization, to 
build up a large enough r es , but slow down late ionization, to de- 



lay the final overlap until z ~ 6 (Wyithe & Loeb 2003; Ciardi et 
al. 2003). Finally, several authors have explored the possibility 
of early partial reionization due to a decaying particle (Chen & 
Kamionkowski 2004; Hansen & Haiman 2004), complemented 
by later full reionization from astrophysical sources. 

The role of small-scale inhomogeneities as sinks of ioniz- 
ing photons has mostly been ignored in this context. Never- 
theless, over a large range of redshifts, the recombination time 
f rec at the mean IGM density is on the order of the correspond- 
ing Hubble time, as illustrated in Figure 1 . Thus the absorption 
of ionizing photons during reionization happens predominantly 
in overdense regions. In hierarchical models like Cold Dark 
Matter (CDM), the smallest structures are the first to collapse 
gravitationally and dominate the photon consumption both dur- 
ing the ionization of a region and afterwards, while balancing 
recombinations. 

When the first sources turned on, they ionized the neutral, 
opaque IGM around them by propagating weak R-type ioniza- 
tion fronts (I-fronts). This type of front moves outward super- 
sonically with respect to both the neutral gas in front of it and 
the ionized gas behind it, so it races ahead of the hydrodynam- 
ical response of the IGM. This process was first described by 
Shapiro (1986) and Shapiro & Giroux (1987), who solved an- 
alytically for the time-varying radius of a spherical I-front sur- 
rounding a point source in the expanding IGM and then used 
this solution to determine when H II regions would grow to the 
point of overlap, thereby completing reionization. In this study 
the effect of density inhomogeneity on the motion of the I-front 
was described by a mean gas clumping factor C = (n 2 )/(n) 2 . 
A clumpy gas has C > 1, which causes the ionized gas to re- 
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FIG. 1 . — Timescales. Hubble time t H (dashed line) and recombination time f re 
the ratio of these timescales (top panel). 
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combine more frequently, increasing the opacity of the H II re- 
gion to ionizing photons, which reduces the flux reaching the 
I-front and slows it down. This approach has formed the basis 
for many more recent semi-analytical treatments of reionization 
(e.g. Haiman & Holder 2003; Wyithe & Loeb 2003; Venkate- 
san, Tumlinson, & Shull 2003). The idea that reionization pro- 
ceeded by the propagation of weak, R-type I-fronts which move 
too fast to be affected by the gas dynamical disturbance they 
create is also the basis for most of the numerical simulations 
of reionization carried out to date (e.g. Razoumov et al. 2002; 
Ciardi, Ferrara, & White, 2003; Sokasian et al. 2004). In par- 
ticular, all numerical studies that add radiative transfer to a pre- 
computed inhomogeneous cosmological density field (i.e. the 
"static" limit) are assuming that there is no significant back re- 
action on the gas 4 . 

The assumption of either a mean clumping factor or the static 
limit to model the effect of density inhomogeneity on cosmo- 
logical I-fronts is not correct even on average, however, unless 
the clumps are either optically thin or absorb only a small frac- 
tion of the ionizing flux. If a clump is self-shielding, then the 
I-front that encounters it will not remain a weak R-type front if 
the size of the clump is larger than its Stromgren length (i.e. the 
length of a column of gas within which the unshielded arrival 
rate of ionizing photons just balances the recombination rate). 
In that case the denser gas of the clump must slow the I-front 
down enough that the disturbed gas inside the clump catches 
up to the I-front and affects its progress. This transforms the I- 
front from supersonic, R-type, to subsonic, D-type and "traps" 
the I-front inside the clump (Shapiro, Iliev & Raga 2004). If 
the clump is gravitationally bound before the arrival of the I- 
front, then the I-front will expel the gas from the clump as a 
supersonic evaporative wind, as long as the clump cannot bind 
photoionized gas with T > 10 4 K. 

The impact of small-scale inhomogeneities on the global I- 
fronts that reionized the universe depended upon the relative 



importance of unshielded and shielded overdense regions and 
their sizes, densities, and abundances. These IGM inhomo- 
geneities can be divided into two major types, both of which 
have been modeled only crudely in the majority of reioniza- 
tion studies. Pre-virialized objects, such as filaments and still- 
collapsing halos, are usually described in terms of the mean 
clumping factor described above. Current semi-analytical mod- 
els of reionization either assume a constant clumping factor 
(Cen 2003; Haiman & Holder 2003; Tumlinson, Venkatesan 
& Shull 2004), a clumping factor derived from linear theory 
(Miralda-Escude, Haehnelt & Rees 2000; Chiu, Fan & Ostriker 
2003; Wyithe & Loeb 2003), or ignore clumping altogether (i.e. 
assume C = 1) (Onken & Miralda-Escude 2004). In practice all 
these approaches are over-simplified since the clumping factor 
of the IGM gas is dominated by the highly-overdense nonlinear 
regions and evolves strongly with redshift. 

Modeling of virialized inhomogeneities in previous studies 
has been even more approximate. An important dividing line 
that separates two distinct populations of virialized halos is that 
defined by the virial temperature, r v j r = 10 4 K. In order for stars 
to form inside halos, the gas must cool below the virial temper- 
ature to become self-gravitating and gravitationally unstable. 
Radiative cooling in a purely atomic gas of primordial com- 
position is ineffective below 10 4 K, however, so "minihalos"- 
halos in the mass range 10 4 M Q <; M <; 10 8 M Q , with virial tem- 
peratures below 10 4 K - are only able to form stars by forming 
H2 molecules, which have the potential to cool the gas below 
the virial temperature, by rotational-vibrational line excitations. 
The H2 that forms in minihalos, though, is easily dissociated by 
UV photons in the Lyman- Werner bands between 11.2 and 13.6 
eV, which are produced in abundance by the first stars, long be- 
fore the ionizing background from such stars is able to reionize 
a significant fraction of the universe (e.g. Haiman, Rees, & 
Loeb 1997; Haiman, Abel & Rees 2000; Ciardi et al. 2000). 
Thus a generic prediction of current structure formation mod- 
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els is a large population of minihalos that are unable to cool and 
form stars. 

In that case, the dominant source of photons for reionization 
would have been the more massive halos (i.e. M J> 10 8 M©) with 
> 10 4 K, in which atomic line cooling is efficient enough to 
enable stars to form. From the point of view of such a source 
halo, all lines of sight will intersect a minihalo at a distance 
less than the mean spacing between sources (Haiman, Abel 
& Madau 2001; Shapiro 2001; Shapiro, Iliev, Raga & Martel 
2003; Shapiro, Iliev & Raga 2004). Thus, the intergalactic I- 
fronts must have found their paths blocked by minihalos in ev- 
ery direction, which trapped the fronts until the minihalos were 
evaporated. 

In Shapiro, Iliev & Raga (2004) and Iliev, Shapiro & Raga 
(2004), we used high-resolution numerical gas dynamical simu- 
lations with radiative transfer to study the encounter between an 
intergalactic I-front and a minihalo in detail, for a wide range of 
conditions expected during reionization. These results yielded 
the number of ionizing photons absorbed per minihalo atom 
during the time between the arrival of the I-front and the evap- 
oration of the minihalo gas, £, as a function of the minihalo 
mass, source flux level and spectrum, and the redshift of the en- 
counter. This is a fundamental ingredient we will need here to 
determine how the presence of minihalos affected global reion- 
ization. 

This trapping of intergalactic I-fronts by minihalos, com- 
bined with the increased recombination rate inside already ion- 
ized regions due to small-scale clumping outside the minihalos, 
may help to explain how reionization could have started early 
and ended late. As the global I-fronts advanced into fresh neu- 
tral regions, they generally encountered minihalos that formed 
at the unfiltered (i.e. not affected by any radiation feedback) 
rate of the universe without reionization. The mass fraction 
collapsed into minihalos in such regions grew over time, from 
8% to 24% to 31% from z = 15 to 9 to 6, so the average number 
of extra photons consumed per atom by photoevaporation must 
also have increased with time. This may have enabled miniha- 
los to slow the advance of the global I-fronts, with increasing 
effect toward late times. Reionization simulations by Ciardi et 
al. (2003) for example, which neglected minihalos, found that 
if one assumes a high escape fraction of ionizing photons from 
the source halos, then the large value of electron scattering op- 
tical depth, T es observed by WMAP can be achieved by the first 
stars in galaxies with mass M £ 1O 9 M as sources. However, 
in this case, reionization is completed far too early. Minihalos 
may have the potential to reconcile this discrepancy, increasing 
the duration of the epoch of reionization and allowing for a sim- 
ilar high value of T es , while postponing the redshift of overlap. 

In this paper, we consider the impact of both minihalos and 
more general IGM clumping in detail, and attempt to quantify 
their effects on the duration of the reionization epoch. Driven 
by measurements of the cosmic microwave background, the 
number abundance of galaxy clusters, and high redshift super- 
nova distance estimates (e.g. Spergel et al. 2003; Eke, Cole & 
Frenk 1996; Perlmutter et al. 1999) we focus our attention on 
the ACDM cosmological model with parameters h = 0.7, f^o = 
0.3, fl A = 0.7, n h = 0.05, cr 8 = 0.87, and n p = 1, where fi 0) H A , 
and Q,b are the total matter, vacuum, and baryonic densities in 
units of the critical density (p C ritX &\ is the variance of linear 
fluctuations filtered on the 8/z _1 Mpc scale, and n p is the index 
of the primordial power spectrum. The Eisenstein & Hu (1999) 
transfer function is used throughout. 



The structure of this work is as follows. In §2 we general- 
ize the approach of Shapiro & Giroux (1987) to account for the 
effect of minihalo evaporation on the time-varying radius of a 
spherical I-front. This will require us to calculate the statisti- 
cally biased abundance of minihalos at the location of the front 
and incorporate the simulation results for the ionizing photon 
consumption rates per minihalo. In §3 we model the global 
progress of reionization by summing the results from Section 
2 over a statistical distribution of source halos, leading to the 
eventual overlap of neighboring H II regions and the comple- 
tion of reionization. Our results and conclusions are given in 
§4. 

2. THE PROPAGATION OF A COSMOLOGICAL IONIZATION 
FRONT ABOUT A SINGLE SOURCE 

2. 1 . Cosmological Ionization Fronts in a Clumpy IGM 

When a source of ionizing radiation turns on in the expand- 
ing, neutral IGM, a weak, R-type I-front propagates outward. If 
the IGM were static, this front would decelerate continuously 
from the moment of turn on, until, within a time comparable to 
the recombination time, it almost reached the size of the Strom- 
gren sphere. This Stromgren sphere is just large enough that the 
total recombination rate of ionized atoms inside it equals the 
ionizing photon luminosity of the central source. At this point 
the I-front drops to the R-critical speed of twice the sound speed 
of the ionized gas, and the front transforms from R-type to D- 
type, preceded by a shock. Thereafter the I-front is affected by 
the dynamical response of the IGM. 

This is not the case, however, in the expanding, average 
IGM. Shapiro & Giroux (1987) showed that, while it is for- 
mally possible to define an "instantaneous" Stromgren radius 
(which grows in time in proportion to the cosmic scale factor), 
the actual I-front generally does not reach this radius. Instead, 
the I-front remains a weak R-type front as long as the source 
continues to shine, and it would not be correct, therefore, to 
describe the cosmological H II region as a Stromgren sphere, 
a misnomer which unfortunately appears in the literature of 
reionization. 

We shall follow the approach of Shapiro & Giroux (1987), in 
which the H II region is bounded by an I-front whose speed is 
determined by the I-front continuity jump condition, which bal- 
ances the outward flux of ionizing photons against the inward 
flux of newly created ions. The flux that reaches the front will 
be determined by solving the equation of transfer between the 
source and the front. We will assume that the IGM is spheri- 
cally symmetric outside the source. For simplicity the I-fronts 
are taken to be "sharp", i.e. the width of the transition between 
the ionized region inside and the neutral region outside the front 
is small compared to its radius. The actual width is compara- 
ble to the absorption mean free path on the neutral side. This 
assumption of small mean free path is generally a good ap- 
proximation for a "soft" Population II (Pop. II) stellar spec- 
trum where most ionizing photons have energies near the ion- 
ization threshold of hydrogen, for which the absorption cross- 
section due to neutral hydrogen is large. However, I-fronts are 
somewhat wider for "hard" spectra like those expected for mas- 
sive Population III (Pop. Ill) stars and the power-law spectra of 
QSOs. In these cases a larger fraction of the ionizing photons 
are at higher energies, corresponding to lower ionization cross- 
sections of neutral hydrogen and helium, and thus our approxi- 
mation of sharp I-fronts is less accurate. 

Adopting this picture, we consider an ionizing source emit- 
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ting N 7 ionizing photons per unit time. We define the comoving 
radius of the ionized region as r/(f) and its comoving volume as 
Vj = 4-nr 3 /3. We are interested in H II regions that at all times 
are much smaller than the scale of the current horizon. The 
jump condition across the I-front is given by a balance of the 
flux of neutral atoms and photoionizing photons. In the frame 
of the front, it can be written as 

n H ,im=fi l F, (1) 

where is the undisturbed hydrogen number density (in 
proper coordinates) on the neutral side of the front, u\ = 
a(drj/dt) is the I-front peculiar velocity, a = 1/(1+ z) is the 
scale factor, and /3, is the number of ionizing photons absorbed 
to create each ionized H atom that emerges on the ionized side. 
In the absence of minihalos, /?, = XeS = 1 + pA(He), which cor- 
rects for the presence of helium, with p = 0, 1 or 2 if He is 
mostly neutral, singly ionized or doubly ionized after the pas- 
sage through the front, and A(He) = 0.08 is the He abundance 
by number with respect to hydrogen. Finally, F is the number 
flux of ionizing photons at the current position of the I-front, 

where S(r,t) is the number of photons emitted by the central 
source which pass through a sphere of comoving radius r per 
unit time, given at r = r t as 

47T 



S{r,, t)=N 1 - —r]a-"(n H ) 2 Ca BX efi 



(3) 



i.e. the number of photons emitted by the source per unit time 
minus the number of recombinations in the current H II re- 
gion volume. Here C is the volume-averaged clumping factor, 
a B = 2.6 x 10~ 13 cm 3 s" 1 is the case B recombination coefficient 
for hydrogen at 10 4 K and is the comoving number density 
of hydrogen in present units, 1.87 x 10~ 7 (fifc/z 2 /0.022) cm -3 . 
In all cases, it is safe to assume that the H II regions are cos- 
mologically small, and hence no ionizing photons are lost to 
redshifting below the hydrogen ionization threshold. 

Combining equations (1) - (3), the evolution of the comoving 
volume of the ionized region V/ is given by 

dVi o dri 1 , n 

= 4 ^7-7T = " ^V 7 -a B C(l + z) 3 4v,. (4) 



dt 



Defining 



dt Xeff«H 



Xeffa B C(«ti) 
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we can write equation (4) in dimensionless form 



dy 
dx 



= \-y{\ + zY 



(5) 



(6) 



where y = V//V s ,i = (n/rs.i) 3 , x = t/t\ and h = l/(a B Cn^) is 
the recombination time of the mean IGM at present (Shapiro & 
Giroux 1987). 

If we define dr = dx/a 3 = (1 +z) 3 dx, equation (6) becomes 
dy =(l+z)- 3 -y, (7) 



dr 

for which a formal solution is 



y{t)-- 



-T(t) 



r(0 



dr'- 



, T(tl) [i+z(t')F 

where f, is the time of source turn-on and for the flat, ACDM 
model 

dr (1 + z) 2 



(8) 



dz Hofi[tto(l+z) 3 + 0\] 1/2 ' 



(9) 



Equation (9) has a solution 

T(z) = K {l-[f!o(l+z) 3 + ft A ] 1/2 }. 



(10) 



where k = 3ff 2 H Qo and the arbitrary constant of integration is 
chosen so that r(z = 0) = 0. At high redshift, before and during 
reionization, we have >> Oa/(1 + z) 3 and the solution (10) 
becomes 

T{z) = n[\-nf(\+z) 312 ]. (11) 
In this limit equation (8) simplifies to 

rrO) 



Equation (12) has an exact analytical solution given by 



n 



(i+zd 3 



tj I 



(12) 

(13) 
dt is 



where rj = 2(1 +z,) 3 / 2 /(3i/ ?i^ /2 ) and Ei(2,x) = f° 
the Exponential integral of second order. The solution in equa- 
tion (13) reduces to the one in equation (10a) of Shapiro & 
Giroux ( 1 987) for a flat, matter-dominated universe with Slo = 1 , 
t = 2/(3//), and thus -q = (1 +z i ) 3 t i /t l . 

2.2. The Average Effect of Minihalo Evaporation on 
Cosmological Ionization Front Propagation 

Having outlined a formalism to describe the expansion of 
ionization fronts in a ACDM cosmology, we next address the 
question of absorption by minihalos. As described in §1, when 
an intergalactic I-front encounters an individual minihalo, it is 
trapped until the minihalo gas is evaporated. For every mini- 
halo atom, this process consumes £ ionizing photons. Suppose 
we consider the average effect of this process on the global I- 
front which moves through a medium comprised of minihalos 
embedded in the IGM. The average speed throughout this com- 
pound medium will be given by a modified I-front continuity 
jump condition which takes account of the additional photon 
consumption due to minihalos. In particular, the quantity /?,■ in 
eq. (1) should now be replaced by the following 

(14) 

where / C oii is the total collapsed baryon fraction (i.e. over all 
halo masses) and / co u,mh is the collapsed fraction of just the 
minihalos. Finally, if £ is the the number of ionizing photons 
consumed per minihalo atom in the encounter between the in- 
tergalactic I-front and an individual minihalo, then £ is the ap- 
propriate average over the distribution of minihalos at the in- 
stantaneous location of the global I-front. Inserting eq. (14) 
into eq. (1) then yields 

F 

n H \u\ = =, (15) 

(1 -/coll)Xeff + [1 +A(He)]/ coU ,MH£ 

where refers to the total H atom density, including both the 
IGM and all halos. 

As usual the flux, F, in this I-front jump condition is deter- 
mined by integrating the equation of transfer over the ionized 
region between the source halo and the I-front. By definition, 
the minihalos originally inside this region do not affect this 
integration, however, since they will already have been evap- 
orated by the passage of the global I-front, thereby returning 
their atoms to the IGM inside the H II region. We assume, for 
simplicity, that the evaporated minihalo gas shares the mean 



A = (1 -/coii)Xeff+ [1 +A(He)]/ co u, MH £, 
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clumping factor the IGM into which it is mixed. In that case, 
the flux at the I-front is given by 

N„-a B C(\+zf (n H ) 2 Xet f (V, - V ) 



F = 



Ana 2 r 



(16) 



where we have been careful now to start our integration of the 
transfer equation from the Lagrangian volume of the source 
halo (i.e. a volume which, when multiplied by the mean den- 
sity, gives the mass of the source halo). 

Combining equations (15) and (16), the evolution of the co- 
moving volume of the H II region in the presence of minihalos 
is then given by 

dVi _ 1 

~dt ~ [(1 -/coii)Xeff + [1 +A(He)]/ coll! M H e] 



-a B C(l+zfn° HXe ff(Vi-V Q ) 



(17) 



where initially V/ = Vo. We note that the solutions in §2.1 in 
the absence of minihalos, including the exact analytical solu- 
tion for Vj(t) and r/(f) in equation (13), will be valid here in 
the presence of minihalos as well, if /3, in equation (14) and the 
clumping factor C are constants, and t\ is redefined as 

tl =(aBCn B ^y. (18) 

In general, both $ and C will not be constant, but nevertheless 
this limit provides a useful check and some insight, as well shall 
see below. 

We adopt a simple model to account for infall when comput- 
ing the flux. The radial coordinate, r If adopted in our formalism 
above is essentially a Lagrangian one, in which we have as- 
sumed that the local mean density of the gas around the source 
is equal to the cosmic mean IGM density. In reality, all mas- 
sive sources are found in overdense regions, due to the gravita- 
tional influence on the surrounding gas. We estimate the map 
between the Lagrangian and Eulerian comoving radii using a 
simple "top-hat" picture, which is a simplified version of the 
model described in Barkana (2004). 

We compute the cross-correlation between a sphere of La- 
grangian radius rj and spherical perturbation of mass equal to 
the source halo mass M s = (47r/3)rg£!oPcrit as 

1 r°° 

a 2 (r ,r,)=— k 2 dkP(k)W(kro)W(kn), (19) 

where P(k) is the initial matter power spectrum, linearly 
extrapolated to the present, and W(x) is the spherical top- 
hat window function, defined in Fourier space as W(x) = 
3 [^j^- -cos(x)x 2 ] , where ro is the Lagrangian radius of the 
source itself. Defining er 2 (ro) = cr 2 (ro,ro), the expected value 
of the linear overdensity of a sphere of Lagrangian radius rj 
about the source is then given by 

1.69D(z)<7 2 (ro,r/) 



+ 1. 



(20) 



D( Zs ) a 2 (r ) 

Here we use Sl to denote the average density within the sphere, 
D(z) is the linear growth factor, z s is the collapse redshift for 
the halo, and the appears since we are defining the over- 
density as p/p instead of p/p - I. For simplicity, we assume 
z = z s , i.e. that the growth of the mass after it collapses due to 
secondary infall, is small. _ 

The linear overdensity Sl can be related to the correspond- 
ing nonlinear overdensity, S, by the standard spherical collapse 



model (e.g. Peebles 1980). In this case, both quantities can be 
expressed parametrically, in terms of a "collapse parameter" 9 

as 



and 



^_ 9 (<9-sin(9) 2 
~ 2(l-cos6») 3: 

2/3 



Sl = 



(6»-sin6») 2/3 +l. 



(21) 



(22) 



These equations define the relationship between the Eulerian 
and Lagrangian comoving radii as 

r,,E = r I 6(r h r Q y 1 / 3 , (23) 

since S given by equation (21) above is the average overdensity 
within an Eulerian sphere of Lagrangian radius rj. 
In the presence of infall, equation (17) becomes 

dVi 1 

dt 



[(1 -/coii)Xeff+ [1 +A(He)]/ coll , MH £] 

Vle 



-a fi C(l+z) 3 4xeff 



(24) 



where now 6 is not the average 5 within the sphere, but rather 
5 at the boundary, and V/ £ = V/<5 _1 is the Eulerian volume. We 
compute 5 as follows. The comoving volume satisfies: 

AV^SiV^ + V,^ S(V,,e) = [V It E + AV,, E ] 8(V ItE + AV ItE ), (25) 

where AVi iE is a small change in the size of the radius. Working 
to first order in AVj jE , this gives 

S = S+V IiE —— . (26) 
dVifi 

We can therefore rewrite equation (24) using the Lagrangian 
volume as 



dV, 
dt 



1 



(27) 



[(1 -/cou)Xeff+ [1 +A(He)]/ coll , MH £] 

^-a B c(i+z) 3 4xef f / ' dv;s(v;r l s(v;f 
A Jv B 

The relevant overdensity that appears in equation (27) is 

i r v > 

<5cium P (V/) = tt^t / dV;S(V;). (28) 

In Figure 2, we show 6, 8, and S c \ ump around peaks of mass 
10 8 M Q and 10 n M Q . Note that 4ium P exceeds 5 over a range 
of radii, because 4ium P is an average in Lagrangian coordinates 
and 5 is an average in Eulerian coordinates. Finally, the flux in 
equation (16) is corrected by a similar factor of <5 c i U mp, yielding 

A>, " Oi B C ( 1 + Z) 3 («£) 2 Xeff (Vl - Vo)4lump(V/) 



F = 



4na 2 i 



(29) 



l.E 



The average number of ionizing photons absorbed per mini- 
halo atom, £, in the process of evaporating all the minihalos at 
the current location of the I-front, r/(f), must now be specified. 
Iliev, Shapiro & Raga (2004) have shown that the number of 
photons per atom absorbed by a minihalo of mass M 7 (in units 
of 10 7 Mq), overtaken by an intergalactic I-front at a redshift 
of z which is driven by an external source of flux Fq, the flux 
in units of that from a source emitting A^ p h = 10 56 s _1 ionizing 
photons per second at a proper distance d of 1 Mpc, i.e. 



d 2 

"Mpc 



8.356 x lOV^m- 2 ' 



(30) 
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FIG. 2. — Overdensites around a source of a given mass, as labeled, versus Lagrangian distance from the center of the source (in units of the Lagrangian radius 
of the halo): the mean S is solid (colored blue in electronic edition), the density at the boundary, <5, is short-dashed (colored red in electronic edition), and <5 c i„ m p is 
long-dashed (colored green in electronic edition). (See the electronic edition of the Journal for the color version of this figure.) 



during its photoevaporation is given by 



(31) 



where <j>\(M) =A(M 7 

7 D+Eiog w F 



B+C\og m M 7 



), <fe(z)= [G + H(l+z)/10] and 



h(Fo) = /^""-'"sio'o Here the factors A-G are dependent on 
the spectrum of the ionizing sources: A = (4.4, 4.0, 2.4), B 
=(0.334, 0.364, 0.338), C = (0.023, 0.033, 0.032), D =(0.199, 

0. 240, 0.219), E = (-0.042,-0.021,-0.036), G=(0.,0.,0.1), H=(l, 

1, 0.9), for the cases in which the ionizing spectrum is taken to 
be a 5 x IQ 4 K blackbody representing Pop. II stars , a QSO-like 
power-law spectrum with slope of -1.8, or a 10 5 K blackbody 
representing Pop. Ill stars, respectively. 

In order to use this result in equation (27), the quantity 
£(M,z,Fq) must first be averaged over the mass function of 
minihalos at rj(t) on the undisturbed side of the I-front as the 
H II region evolves with time. Since £(M,z,Fq) depends on 
Fo(t), which also depends on r/(f) according to equation (29), 
equations (27) and (29) are coupled and must be solved simulta- 
neously. We shall consider three analytical approximations for 
the minihalo mass function when averaging £(M,z,Fq) to ob- 
tain £. The first two approaches, described in § 2.2. 1, are based 
on the well-known Press-Schechter (PS) approximation for the 
mass function averaged over all space at a given redshift (Press 
& Schechter 1974). We shall refer to these, which depend upon 
Z and Fq, but not upon rj(t ) or the source halo properties, as "un- 
biased minihalo" averages. The third approximation, described 
in § 2.2.2, is based upon an extension of the PS approach which 
takes account of the spatial correlation between the minihalos 
and the central source halos, as described by Scannapieco & 
Barkana (2002). In this last approximation, which we refer to 
as the "biased minihalo" average, £ not only depends upon z 
and Fq, but also on r/(f) and the source halo mass. 



2.2.1. The Average Photon Consumption Rate for Unbiased 

Minihalos. 



We begin by defining the average photon consumption rate 
per minihalo atom 



£,nb,\(z,F Q ) 



dn(M.z) 
dM 



M£(M,z,F ) 



f M ~ dM dn{M ' z) M 



(32) 



where ' fa ^r ) is the PS mass function of halos, if we assume 

dM ' 

that the minihalos at a given redshift z just formed at that red- 
shift. If, on the other hand, we assume that minihalos at z had a 
distribution of formation redshifts zj, with z/ > Z, then 



iz: dM ir^f d -^M^M,z f ,F 0) 



£,nb,i(z,F Q ) : 



iMi dM ir dz f 



d 2 n(M,z { ) 



M 



(33) 



In both these equations, the limit M m m is the minimum mini- 
halo mass (which we assume here to be the Jeans mass at that 
epoch), while M max = M{T VU = 10 4 /T) is the halo mass at that 
epoch for which T v \ r = 10 4 K. In equation (33) we have ap- 
proximately accounted for the distribution of minihalo forma- 
tion times, by taking the derivative of the mass function, which 
glosses over the fact that the change in this function at a given 
mass M includes both a positive contribution from halos whose 
masses have increased to M from lower values, as well as a 
negative contribution from halos whose masses have increased 
from M to higher values. The error introduced by this approx- 
imation is small, however (eg. Kitayama & Suto 1996), and is 
justified given the other uncertainties involved. 

2.2.2. The Average Photon Consumption Rate for Biased 
Minihalos. 

To calculate the biased distribution of minihalos about a 
given source, we employ an analytical formalism that tracks 
the correlated formation of objects. Our approach, described 
in detail in Scannapieco & Barkana (2002), extends the stan- 
dard PS method using a simple approximation to construct the 
bivariate mass function of two perturbations of arbitrary mass 
and collapse redshift, initially separated by a fixed comoving 
distance (see also Porciani et al. 1998). From this function we 
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can construct the number density of minihalos of mass M that 
form at an initial redshift z at a comoving distance r from the 
source halo of mass M s and formation redshift z s - 



dM {M, Z ,r\M s ,z s ) *L {MsjZs) > 



(34) 



dM, 



where j^-(M s ,z s ) is the usual PS mass function and 



d L n 



dMdM ~ (M,z,M s ,z s ,r) is the bivariate mass function that gives 
the product of the differential number densities at two points 
separated by an initial comoving distance r, at any two masses 
and redshifts. Note that this expression interpolates smoothly 
between all standard analytical limits: reducing, for example, 
to the standard halo bias expression described by Mo & White 
(1996) in the limit of equal mass halos at the same redshift, and 
reproducing the Lacey & Cole (1993) progenitor distribution in 
the limit of different-mass halos at the same position at different 
redshifts. Further detaled checks of the method against N-body 
simulations were presented in Iliev et al. (2003) and Scanna- 
pieco & Thacker (2005). Note also that in adopting this defini- 
tion we are effectively working in Lagrangian space, such that 
r is the initial comoving distance between the perturbations. As 
a shorthand we define ^ £ (z,r) = ^(M,z,r\M s ,z s ). With this 
definition the biased values of the average photon consumption 
per minihalo atom corresponding to the two unbiased cases in 
equations (32) and (33), respectively, are 



£, h . l (z,F Q ,ri,M s ) = 



JMldM a -^(z, ri )MaM,z,F ) 



, dn mh 



L 



■Win 



and 



dU 

dM d -^r(z,rj)M 



_ iMldMirdzf^izf^M^M^Fo) 



(35) 



£,b,2{z,Fo,r h M s ) = 



lMldMf-dz f ^(z f , ri )M 



(36) 

with n the Lagrangian radius of the I-front. 

We can anticipate how important this bias effect is likely to 
be in determining the average minihalo consumption rate by 
considering its impact on the minihalo collapsed fraction in the 
vicinity of a given source halo. In Figure 3, we plot the av- 
erage (unbiased) minihalo collapsed fraction versus the biased 
value, which varies with distance from the source halo, for two 
halo masses 1O 8 M and 1O H M , at redshifts z = 7 and z = 15. 
These halo masses and redshifts illustrate the range of behavior 
expected during reionization. In terms of the Gaussian-random- 
noise initial conditions for our ACDM cosmological model, 
10 8 M Q halos correspond to fluctuations that are 3.2 a (1 .6 a) at 
z = 15 (7), respectively while 1O U M halos are 6.4 a (3.2 a) at 
these same redshifts. 

The distance between sources and minihalos is measured in 
terms of the comoving Lagrangian (i.e. unperturbed) radius of 
a given mass shell surrounding the central source, in units of ro, 
the Lagrangian radius of the source. According to this figure, 
the bias can be significant for minihalos located in the range 
1 <J r/ro <J 10. For typical source halos, in fact, the biased col- 
lapsed fraction in the neighborhood of the source hardly de- 
clines with increasing redshift, in contrast with the unbiased 
collapsed fraction which declines by a factor of more than 3 
between z = 7 and z = 15. 

2.3. The Ionizing Photon Luminosity of the Central Source 



The total ionizing photon output of a source, N y , and its time 
evolution depend on the mass of the host halo M s , photon pro- 
duction per stellar baryon A 7 ,-, star-formation efficiency and 
ionizing photon escape fraction / esc . We then define the total 
ionizing photon output per source atom that escapes the source 
halo as: 

/ 7 =/*/esc^, (37) 

thus the source emits a total of N 7 = f^MVLb / (pm p ) during its 
lifetime, where /j,m p is the mean mass per atom. 

The star-formation efficiency /* and ionizing photon escape 
fraction / esc are highly uncertain in general, and even more so 
for the high-redshift galaxies responsible for reionization. Their 
estimated values vary by several orders of magnitude between 
different observational and theoretical estimates (Leitherer et 
al. 1995; Ricotti & Shull 2000; Heckman et al. 2001; Stei- 
del et al. 2001; Tan & McKee 2001). For simplicity, and 
since the principle aim of this investigation is to show the ef- 
fect of small-scale structure rather than model the sources in 
detail, we assume that each source produces a fixed number 
/ 7 of ionizing photons which escape from the source galaxy 
per atom in the source during the source's lifetime. The ion- 
izing photon production per atom for Pop. II low-metalicity 
stars with a Salpeter IMF is N t = 3000- 10000 (Leitherer et 
al. 1999). Zero-metalicity, massive Pop. Ill stars, on the other 
hand, are estimated to produce values of Nj that rise sharply 
with mass from 25,000 to 80,000 as stellar mass increases from 
10 M to 50 M Q , then gradually reach a peak of 90,000 at 120 
M Q , and finally decline slowly to 80,000 by 500 M Q (Schaerer 
2002; Tumlinson, Venkatesan & Shull 2004). While individ- 
ual Pop. Ill stars have higher values of / 7 than Pop. II stars of 
the same mass, a non-trivial part of the increase from Pop. II to 
Pop. Ill quoted above reflects the fact that the assumed Pop. II 
IMF has many low-mass stars, which are inefficient ionizing 
sources, while the Pop. Ill IMF is often hypothesized to con- 
tain only massive stars. 

Assuming, conservatively, that Nj > 4000 for Pop. II stars 
and A 7 ,- > 25,000 for Pop. Ill stars, and taking moderate 
fiducial values for the photon escape fraction and star for- 
mation efficiency of / esc =0.1 and /* = 0.1, yields / 7 > 
(40,250)(/ esc /0.1)(/*/0.1) for (Pop. II, Pop. Ill), respectively. 
We shall further assume that the time-dependence of this ioniz- 
ing photon output is characteristic of a starburst with a photon 
luminosity 

■ 1 Mfl h ( 1 * < 



AW- 



a - 



x 



ay 



t >t s 



(38) 



a jjbm.pt s 

where a = 4.5, i.e. we assume that the source is steady for a 
time t s , after which the photon flux decreases as power of time 
(Haiman & Holder 2003). Here t s is the characteristic time for 
a source to fade, essentially the typical source lifetime. For in- 
dividual massive stars t s « 3 Myr, but starbursts could in prin- 
ciple last significantly longer, thus we consider both the cases 
f s = 3 Myr and t s =100 Myr. We assume that the He correc- 
tion Xeff is L08 for the softer Pop. II spectrum and 1.16 for the 
harder Pop. Ill and QSO spectra. 

2.4. Results for Individual H II Regions 

We present the results of our numerical solution of the spher- 
ical I-front evolution equations from §2.2 for two illustrative 
cases. Sample results for two source halos, a 1O 8 M halo that 
turns on at z = 15 and a 1O U M halo that turns on at z = 7, 
both with / 7 = 250, t s = 3 Myr, and Pop. II spectra, are given 
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FIG. 3. — Biased collapsed fraction of baryons in minihalos, / co ll,MH, as a function of the Lagrangian distance from the source halo (in units of the source halo 
Lagrangian radius) for sources of masses 10 8 Mq and 10"Mq and redshifts z = 15 and z = 7, as indicated (solid). For reference we also show the unbiased collapsed 
fraction of baryons in minihalos at the corresponding redshifts (dashed). (See the electronic edition of the Journal for the color version of this figure.) 
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FIG. 4. — Bottom: The evolution of the Lagrangian volume of the H II region about a single source of (left) mass 10 Mq that turns on at zt = 15, or (right) mass 
10" A/0 that turns on at z = 7. Both sources have Pop. II stellar spectra and lifetimes of t s = 3 Myr, during which they produce a total of / 7 = 250 photons/atom. 
Shown are the cases of no minihalos (solid), unbiased minihalos (dotted), and biased minihalos (dashed) for IGM clumping factors (top to bottom in each case) 
C = (i.e. no recombinations in IGM gas), 1 (mean IGM), and 10 (clumped IGM). V max is the maximum ionized volume reached during the lifetime of the source in 
the C = case with no minihalos, as defined in the text. Top: Ratios of the ionized volumes with unbiased and biased minihalos to the no minihalo case, as labeled, 
for C = 1 (solid) and 10 (dashed). (See the electronic edition of the Journal for the color version of this figure.) 



9 




FIG . 5 . — Evolution of individual H II regions with the same source parameters as in Figure 4. Top: The correction factor ft/Xeff due to minihalos for the number 
of ionized photons consumed per atom that crosses the I-front, for biased (dashed) and unbiased minihalos (dotted) for C = 0, 1 and 10. Bottom: Comoving radius of 
the H II region for no minihalos (solid), unbiased (dotted) and biased minihalos (dashed) for C = 0, 1 and 10. (See the electronic edition of the Journal for the color 
version of this figure.) 
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FIG. 6. — Evolution of the dimensionless ionizing photon flux Fo(') at the current position of the I-front. Same notation as in Fig. 4. Bottom set of (initially- 
overlapping) curves are for C = 0, middle set - for C = 1 and top set - for C = 10, respectively.(See the electronic edition of the Journal for the color version of this 
figure.) 
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FIG. 7. — Evolution of the IGM clumping factor in ACDM from numerical N-body simulations, for the gas outside halos. 
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shift based on the PS formalism. We then calculate the evolu- 
tion of the H II region created by each source halo, as discussed 
in section §2.2 and §2.3. Finally, we add the volumes of all 
these H II regions. This gives the total ionized mass fraction at 
each redshift, according to 



dz! 



dMdz! 



-V,(z,z',M,N^, (39) 



in Figures 4-6. In both cases we display results for a variety of 
clumping factors (C = 0, 1, 10) and successive approximations 
of no minihalos, unbiased minihalos [as given by eq. (32)], and 
biased minihalos [as given by eq. (35)]. In Figure 4 we show 
the evolution of the ionized volume. All volumes are normal- 
ized to V max = N-y/rfifj, i.e. the volume ionized if there were no 
minihalos and no recombinations in the gas. 

The increase of the clumping factor from C=ltoC=10 
significantly decreases the maximum ionized volume achieved 
by the H II regions, with and without minihalos, especially 
at higher redshift, as shown for the z= 15 case in Figure 4. 
For similar reasons, H II regions around source halos in over- 
dense areas of the IGM must also be smaller than those in 
the mean IGM. For C = and no minihalos, Vmax/Vb = / 7 , so 

1 /3 

>wAo = f 7 ■ Since more realistic cases, with C> 1 and 3 L Time-Dependent Clumping Factor of IGM Outside Halos 
minihalos, all have V/ < V max , it must be true that rj/ro < f^ 3 



where d ^j^f } is the PS distribution of the source halos, and 
we assume here that source halos have masses M s > M s , m ; n = 
M(10 4 K), the mass of halos with r v ; r = 10 4 K. The universe is 
fully ionized when the H II regions overlap, which corresponds 
to//,M = 1. 
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in general. According to Figure 2, the overdensity <5 c i um p > 1 
outside the source halo for all r/ro<> 20, while <5 c i ump > 2 for all 
r/r„ <J 4 (6), for source halos of mass 10 11 (10 8 ) M Q . As such, 
the IGM recombination correction to the ionized volume at any 
epoch must be significantly enhanced by this local overdensity 
for any / 7 <; 1000. In short, for a realistic range of f 1 values, 
cosmological H II regions from stellar sources are generally not 
larger than the infall regions associated with their source halos. 
This is true with and without minihalos. 

Next we consider the effect of adding minihalos. According 
to Figure 4, for the 10 8 M Q source halo, the unbiased minihalo 
distribution decreases the ionized volume by ~ 20% compared 
to the no minihalos case, while when we account for the mini- 
halo bias about the source, the ionized volume is decreased by 
a factor of 2 relative to the no minihalos case. For the 1O H M 
source halo, the net ionized volume when minihalos are present 
(biased or not) is about 65% of the volume in the case without 
minihalos. 

In the top panel of Figure 5, we plot the factor /3,/Xeff by 
which minihalo evaporation boosts the number of ionizing pho- 
tons consumed at the I-front per atom that crosses the front (in 
the IGM and in minihalos combined). In the bottom panel of 
this figure, we plot the comoving (Eulerian) radius of the cor- 
responding H II regions. For the 10 8 M Q source halo at zi = 15, 
ignoring the minihalo bias (bottom lines) seriously underesti- 
mates the photon consumption by a factor of ^ 2 as compared 
to the biased minihalos (top lines). Furthermore, the overall ef- 
fect of adding minihalos is to increase the photon consumption 
by <~ 100%. For the 10 n MQ halo at z, = 7, there is a similar in- 
crease in the photon consumption, but little difference between 
the biased and unbiased results. 

Finally we plot the dimensionless flux Fo at the current po- 
sition of the I-front, in Figure 6. This is important as a check 
of our assumptions, which incorporate simulation results for 
minihalo evaporation for a range of fluxes, 10~ 2 < Fo < 10 3 . 
The value of Fo has a significant impact on the ionizing photon 
consumption of minihalos, which is higher for higher values 
of Fo. In both cases, the flux starts fairly high (Fo ;> 100), but 
drops to Fo <~ 1 by the time the source starts to fade. At first this 
drop is due mainly to geometric dilution, but later, recombina- 
tions in the ionized volume accelerate this decrease according 
to equation (29). 

3. TOWARDS A MORE GLOBAL PICTURE 

To construct models of the global reionization process, we 
first calculate the number density of source halos at each red- 



As in the individual-source results above, we calculate the 
evolution of the ionized mass fraction for the cases of no 
minihalos, unbiased minihalos, and biased minihalos, and for 
clumping factors of C = 0, 1 and 10. In addition, we consider 
the more realistic case of a clumping factor that evolves in time 
as more and more structure forms, which was obtained from 
numerical N-body simulations by the particle-particle/particle 
mesh (P 3 M) method, with a computational box size of 1 co- 
moving Mpc with 128 3 particles and 256 3 cells, corresponding 
to a particle mass of 2 x 10 4 M Q (see Shapiro 2001 and Iliev 
et al. 2003 for details on the simulations). This box size was 
chosen so as to resolve the scales which contribute most of the 
clumping - on smaller scales the gas would be Jeans smoothed, 
while on larger scales the density fluctuations are still linear and 
do not contribute much to the overall clumping factor. The re- 
sult for the IGM clumping factor is plotted in Figure 7. This 
clumping factor excludes the matter in collapsed halos, since as 
we discussed above, these are self-shielded and we treat them 
separately. The evolution of this IGM clumping factor with red- 
shift is well-fit by 



C(z) = 17.6<T 01(te+00011 < 



(40) 



3.2. The Global Consumption of Ionizing Photons During the 
Epoch of Reionization 

The reionization of the universe was complete when the vol- 
ume of all the H II regions at some epoch equaled the total 
volume. We call that the epoch of overlap, at redshift z ov . Gen- 
eralizing Shapiro & Giroux (1987) we define a useful dimen- 
sionless ratio of total number of ionizing photons emitted per 
hydrogen atom in the universe until overlap at z = z m given by, 

f' ( ^ ,d 2 n l !(M„z')dz' 



C°v - 



f 7 Q b 



n^jimpts 



A/».„ 



Hz) 



(41) 

where d 2 n®(M,z)/dMdz is the comoving differential number 
density of the source halos of mass M formed at redshift z. 

3.3. Electron Scattering Optical Depth Through the Reionized 

Universe 

For any given reionization history, the mean optical depth 
along a line-of-sight between an observer at z = and a red- 
shift z due to Thomson scattering by free electrons in the post- 
recombination universe is given by 

f° / , dt 
T es (z) = ccr T / dz'n e (z')— , (42) 
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where aj = 6.65 x 10~ 25 cm 2 is the Thomson scattering cross- 
section, c is the speed of light, and n e (z) is the mean number 
density of free electrons at redshift z, given by 

« e (z) = «e,0(l+z) 3 //,M(z), (43) 

where //,m(z) is the ionized fraction of the IGM at redshift z, and 
"e,o = "nXeff- For comparison with the value of r es between us 
and the surface of last scattering inferred from measurements 
of the polarization of the CMB, we should evaluate T es (z) at 
Z = Ziec, the redshift of recombination. For this paper we assume 
that He is ionized to He II at the same time as H is ionized. We 
make the reasonable approximation that p = 2 (1) in \eS f° r 
z < 3.5 (z > 3.5). In fact, the reionization of helium to He III 
at z ~ 3 - 4, inferred from observations of quasar absorption 
spectra, has only a small effect on the total electron-scattering 
optical depth and can be ignored for most purposes. In the case 
of fiM = const between us and redshift z (e.g. for an IGM fully 
ionized and free of minihalos since z ov , f^ M = 1 for z < z ov ), the 
integral in equation (42) has a closed analytical form 

2ca T n bPcnt ^ fm | [0q(1 +z)3 + 0a]V 2 _ x } ; (44) 



7"es(z) = ■ 



3H Q Li H m p Q 

where fi H = 1 +4A(He) = 1.32, the mean molecular weight per 
H atom. For f 1M = 1(0) for z < z ov (z > z ov ), r es > 0.04 for 

Zov > 6. 

4. GLOBAL REIONIZATION MODELS: RESULTS AND 
CONCLUSIONS 

4.1. Fiducial Model 

We first consider a fiducial model, Case A, that assumes that 
the sources are Pop. II starbursts, with / 7 = 250 and a lifetime of 
t s = 3 Myr. The evolution of the ionized mass fractions for Case 
A, with and without minihalos and IGM clumping, is shown in 
Figure 8. At z J> 15 there are very few minihalos. Thus in the 
C = case, the no-minihalos model and the model with unbi- 
ased minihalos are almost equivalent. The few minihalos that 
do exist at such early times, however, are highly biased towards 
the sources, giving rise to a substantial difference between the 
biased and unbiased models. Similar trends are visible in the 
C = 1 case with minihalos having an appreciable impact only 
in the biased case. On the other hand, in the C(z) model, re- 
combinations in non-virialized structures slow reionization to 
the point where a significant average minihalo collapse frac- 
tion builds up, giving rise to a shift (i.e. delay) of Az <~ 2 be- 
tween the no minihalos and unbiased minihalo models. In this 
case, the minihalo correction increases with time relative to the 
case without minihalos, as expected since the unbiased / C o11,mh 
grows with time. This has the effect of slowing the rate of in- 
crease of the ionized volume more and more over time. This 
trend would help reconcile an early onset of reionization with a 
late epoch of overlap. When minihalo bias is taken into account, 
however, the minihalo correction is just as large at early times 
as at late times. At late times the biased and unbiased minihalo 
curves almost converge, since the clustering of small objects is 
very weak at these times and bias has a minimal effect. 

The results for z ov , Cov and the net T es , T es (z rec X f° r our fidu- 
cial model, Case A, are shown in Table 1 and in Figure 9 for 
Tes(z). When there are no minihalos and no recombinations in 
the IGM, Zov = 15.2 and T es = 0.186, which requires just one 
photon per atom (i.e. £ ov = 1). The presence of minihalos by 
themselves delays overlap until z ov =14, while increasing the 
global photon consumption by a factor of ~ 2 and decreas- 
ing the optical depth to 0.169. The IGM clumping by itself 



delays overlap until z ov = 8.1(9.5) for C = 10[C(z)], increas- 
ing the global photon consumption to 23(14) and decreasing 
r es to 0.090(0.114). When the effects of both minihalos and 
redshift-dependent IGM clumping are included, overlap is fur- 
ther delayed until z ov = 7.2, in rough consistency with the re- 
sults from the Gunn-Peterson trough observations, while the 
electron-scattering optical depth decreases to 0.089, somewhat 
below the 1 — a WMAP limit, and the global ionizing photon 
consumption rises to ( m = 33. 

4.2. Impact of Model Uncertainties and Discussion 

In order to estimate the effects of small-scale structure under 
different assumptions, we consider several cases in addition to 
our fiducial model (Case A). Case B is the same as our fidu- 
cial model, but assumes that the sources are longer-lived, with 
t s =100 Myr. Case C is also the same as our fiducial model, 
but assumes that the formation times of the minihalos are dis- 
tributed in redshift, by replacing equation (32) with equation 
(33) for unbiased minihalos and equation (35) with equation 
(36) for biased minihalos as discussed in sections § 2.2.1 and 
§ 2.2.2. Cases D, E and F are like Case A, except with f 7 = 40 
(Case D), 150 (Case E) and 500 (Case F). Finally, Case G as- 
sumes Pop. Ill sources with a t s = 3 Myr lifetime and f~ f = 250. 

The values of z ov , Cov and r es for each of these cases are sum- 
marized in Table 1 . In all cases the clumping of the IGM gas a 
strong effect on the duration of reionization by significantly in- 
creasing recombination rates and the ionized photon consump- 
tion. Thus increasing C from 1 (no clumping) to 10 delays the 
overlap by Az « 4-6 and increases the global ionizing photon 
consumption ( OY by a significant factor, between 6 and 17. The 
electron-scattering optical depth r es decreases correspondingly, 
from ~ 0. 1 1 — 0. 17 (consistent with current WMAP limits) for 
C = 1, to - 0.04-0.1 for C = 10. Interestingly, the cases with 
evolving clumping factor C(z) yield epochs of overlap similar 
to those of the C = 10 cases, while at the same time noticeably 
increasing T es , by Ar es ~ 0.02 to ~ 0.09-0.1 1, reaching better 
agreement with the observational limits from WMAP. 

The presence of minihalos additionally delays reionization, 
by up to Az « 2.5. This effect is strongest for short-lived 
sources, and almost disappears for very long-lived sources 
(Case B). This is because in the long-lived case, the same total 
number of photons are produced over a longer time, and thus 
the typical flux responsible for evaporating minihalos is lower, 
leading to more efficient photoevaporation in terms of ioniz- 
ing photon consumption. Similarly, in the short-lived case, the 
photon consumption by minihalos typically raises the global 
photon consumption per atom until overlap, £ ov , by a factor of 
~ 2, while in Case B, the £ ov increases range from 10 to 60%, 
depending on the clumping factor. 

Turning to the question of our assumptions about minihalo 
formation times, we find that the differences between cases A 
and C are negligible in all cases. Thus we can generally assume 
that the minihalos just formed at the redshift of consideration 
as opposed to formation times that are distributed in redshift. 

What is the effect of varying the efficiency for ionizing pho- 
ton release? In Case D, in which reionization is caused by 
Pop. II sources with / 7 = 40, we find that overlap is achieved be- 
fore z = 6 only when there is no gas clumping (C = 1), even if no 
minihalos are present. Hence, / 7 (which depends upon the as- 
sumed photon production per stellar baryon, star-formation effi- 
ciency per halo baryon and the ionizing photon escape fraction) 
should be significantly larger than this assumed value (e.g. by 
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FIG. 8. — Global reionization. Sources with / 7 = 250 and lifetime of t s = 3 Myr (Case A) are assumed, (bottom panel) Decimal logarithm of ionized mass (or 
Lagrangian volume) fractions (i.e. corresponds to overlap) for the cases of no minihalos (solid), unbiased minihalos (dotted), and biased minihalos (dashed) for 
IGM clumping factors (top to bottom in each case) C = 0, 1 , and C(z) (clumped IGM). (upper panel) Ratios of the ionized volume fractions in the presence of biased 
minihalos and with no minihalos for C = 0, 1 and z-dependent.(See the electronic edition of the Journal for the color version of this figure.) 




FIG. 9. — Global reionization. Integrated optical depth due to electron scattering t cs vs. redshift. Same notation as in Fig. 8. Top (dotted) curve shows the optical 
depth produced assuming complete ionization out to the corresponding redshift. Horizontal lines indicate the best-fit and 1 - a uncertainties of the first-year WMAP 
result, t cs = 0.17 ± 0.04. (See the electronic edition of the Journal for the color version of this figure.) 
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Table 1 

Global reionization results for ionizing photon consumption, £ ov overlap epoch, z ov a , and electron scattering 

OPTICAL DEPTH, T. 



No MH 






C = 






C= 1 






C= 10 




— 


C = C(z) 




Case 




— 






Cov 


Zov 


7" 


— 




Tes 




Zov 


Tes 


A 


250 


1 


15.4 


0.186 


2.1 


14.2 


0.168 


23 


8.1 


0.089 


14 


9.5 


0.1 14 


B 


250,long life 


1 


13.6 


0.157 


1.8 


13 


0.144 


13 


9 


0.096 


9 


9.8 


0.108 


C 


250,z-distr. 


1 


15.4 


0.186 


2 


14.2 


0.168 


24.5 


8.1 


0.090 


15 


9.4 


0.114 


D 


40 


1 


11.3 


0.137 


2.7 


9.3 


0.107 






0.045 






0.048 


E 


150 


1 


14.4 


0.173 


2.3 


12.7 


0.151 


24 


6.5 


0.072 


21 


7.0 


0.088 


F 


500 


1 


17 


0.205 


2. 


15.5 


0.192 


20 


10.5 


0.118 


10.0 


12 


0.151 


G 


250 


1 


15.4 


0.186 


2.1 


14.2 


0.168 


23 


8.1 


0.089 


14 


9.5 


0.114 


MH, no bias 




























A 


250 


1.2 


15 


0.183 


2.5 


13.6 


0.164 


39 


6.7 


0.077 


30 


7.5 


0.099 


B 


250 Ion? life 


1.2 


13.5 


0.155 


2 


12.5 


0.142 


14 


9 


0.096 


9 


9.8 


0.107 


C 


250,z-distr. 


1.2 


15 


0.183 


2.6 


13.7 


0.163 


43 


6.6 


0.076 


34 


7 


0.095 


D 


40 


1.3 


11 


0.131 


4.4 


8 


0.095 






0.043 






0.044 


E 


150 


1.2 


14 


0.169 


3 


12 


0.144 


34 


5.5 


0.059 


34 


5.5 


0.072 


F 


500 


1.2 


16.5 


0.203 


2.1 


15.5 


0.189 


30 


9.5 


0.108 


16.0 


11 


0.141 


G 


250 


1.1 


15.2 


0.185 


2.2 


14 


0.166 


31 


7.3 


0.082 


22 


8.3 


0.106 


MH, w/bias 




























A 


250 


1.9 


14 


0.169 


4.3 


12.4 


0.148 


40 


6.6 


0.072 


33 


7.2 


0.089 


B 


250,long life 


1.6 


13 


0.148 


2.3 


12.3 


0.138 


14 


9 


0.095 


10 


9.7 


0.106 


C 


250,z-distr. 


2.1 


14 


0.169 


4.5 


12.4 


0.148 


44 


6.6 


0.072 


35 


7 


0.087 


D 


40 


2.0 


10 


0.117 


5.1 


7.5 


0.085 






0.043 






0.044 


E 


150 


2.0 


13 


0.154 


4.6 


11.0 


0.129 


34 


5.5 


0.057 


34 


5.5 


0.067 


F 


500 


1.8 


15.3 


0.192 


3.4 


14.5 


0.176 


37 


9 


0.100 


22.0 


10 


0.126 


G 


250 


1.5 


14.7 


0.177 


3.2 


13 


0.157 


32 


7.2 


0.079 


24.5 


8 


0.099 



"When no value of Zov is given, overlap was not achieved by z = 5.5. 



replacing the Salpeter IMF by a top-heavy one). For / 7 = 150 
(Case E), the values of z ov and r es without minihalos are essen- 
tially equivalent to those in our fiducial Case A when biased 
minihalos are included, for all values of the clumping factor. 
The quantity £ ov is smaller in Case E than in Case A, however. 

Since Case A has too small a value of r es to satisfy the 
WMAP constraint, we consider a case with higher efficiency 
for photon release / 7 = 500 (Case F), to increase r es by mak- 
ing reionization earlier. With evolving clumping factor C(z) 
and biased minihalos, Case F results in r es = 0.126, which 
is marginally consistent with the \-a WMAP constraint, but 
z ov = 10, too early to match quasar observations. 

Finally, we illustrate the difference between Pop. II and 
Pop. Ill ionizing source spectra in Case G, for which massive 
Pop. Ill stars provide / 7 = 250. The only difference between 
Case G and Case A is the source spectral shape, which results in 
somewhat more efficient evaporation of the minihalos. Hence, 
overlap is slightly earlier, at z ov = 8, and £ ov is a bit smaller, at 
25, while r es is a bit higher, at r es ~ 0.1, so there is still not a 
good agreement with both the quasar and CMB constraints. If 
we reduce / 7 somewhat, for this Pop. Ill case, the overlap red- 
shift would drop, but so would T es , so it would not help to match 
these observations simultaneously. We conclude from this that 
this problem is not very sensitive to our choice of source spec- 
trum. 

In addition, we checked the dependence of our results (for 
our fiducial Case A) on our assumption that the halo mass func- 
tion is well-described by the PS distribution, by replacing it 
with the Sheth and Tormen mass function (e.g. Sheth & Tormen 
2002, hereafter ST) for the case without bias (our bias formal- 
ism does not apply to ST). The result is that noticeably more 
ionizing sources (i.e. more massive halos) exist at early times, 
which makes the electron scattering optical depth slightly larger 



and thus a bit closer to the one-sigma WMAP range. However, 
the effect is rather marginal, and, in the presence of minihalos 
and evolving IGM clumping the overlap occurs at almost ex- 
actly the same time using either PS or ST. In the case with no 
bias, at least, we conclude that using ST instead of PS does not 
make a significant difference for the current reionization con- 
straints. When bias is included we expect that the differences 
will be even smaller, since biased minihalos cluster strongly 
around the sources, thus screening them and consuming more 
ionizing photons than in the unbiased case, which should make 
the higher abundance of sources early on less important. It 
should also be noted that currently ST is not well tested at high 
redshifts, hence it is not clear if it is more precise than the stan- 
dard PS model. In fact, the few simulations of high-z small- 
scale structure formation that exist seem to agree with PS rea- 
sonably well (Shapiro 2001; Jang-Condell & Hernquist 2001; 
Cen et al. 2004). 

There are some additional effects which we have neglected 
here which we intend to treat in a future paper. When the lumi- 
nosity of a central source halo inside an H II region decays, the 
size of the H II region shrinks. Gas that was once ionized and 
inside the H II region before this, can later find itself outside 
the H II region, where it no longer receives ionizing radiation 
from the central source. In our current treatment, such gas is 
assumed to recombine and must be ionized again before reion- 
ization of the universe is complete. Formation of minihalos in 
such gas will be suppressed by the Jeans mass filtering of the 
IGM (Shapiro, Giroux & Babul 1994; Gnedin 2000b) and, to a 
lesser extent, by the relic entropy of the gas after it cools and 
recombines (Oh & Haiman 2003). A similar effect occurs if 
and when additional sources of heating of the IGM are present 
beyond the ionizing sources we have considered, such as an X- 
ray background (Machacek, Bryan & Abel 2003; Oh & Haiman 
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2003; Glover & Brandt 2003). Our neglect of these possibili- 
ties for the partial suppression of the minihalo population may 
overestimate the contribution of minihalos to ionizing photon 
consumption before final overlap. However, this effect is off- 
set by another, perhaps more important, effect which we have 
also neglected. We have considered the clustering of minihalos 
around source halos, but we should also consider the clustering 
of source halos in space and in time. The latter effect makes it 
possible for new sources to ionize the pre-existing ionized vol- 
ume of an H II region created by one source, preventing the H 
II regions from ever shrinking. In fact, if current simulations of 
reionization are correct in this regard, then this effect is likely 
to be important, since simulations find that ionized zones tend 
to grow over time and eventually overlap, rather than shrinking 
and disappearing and then being replaced by independent H II 
regions that form elsewhere (e.g. Ricotti, Gnedin & Shull 2002; 
Ciardi, Ferrara & White 2003; Sokasian et al. 2004). A simi- 
lar point has been made recently by Furlanetto, Zaldarriaga & 
Hernquist (2004). In that case, the suppression of minihalo for- 
mation in pre-ionized zones which recombine and must be ion- 
ized again, mentioned above, is unimportant, since gas which 
is ionized once, remains ionized by the continuous exposure to 
new sources over time. We will treat these effects in a future 
paper. 

It would appear then that our fiducial model Case A with 
minihalos is the best candidate for helping to resolve QSO con- 
straints on z ov with the WMAP measurements of r es . For al- 
though the model still falls short of reproducing the observed 
r es , minihalos nevertheless form on average in small numbers 
at early times and in large numbers at late times, thereby ex- 
tending reionization and accumulating scattering optical depth 
before overlap. However, this behavior changes dramatically 
if one accounts for the strong clustering of high-z minihalos. 
As shown in Figure 3, the minihalo collapsed fraction in the 
neighborhood of the typical source halos is significantly higher 
than the spatially-averaged collapsed fraction of minihalos at 
the same epoch. Hence, although the total number of minihalos 
increases with time, the number around sources remains largely 
fixed, pushing the entire process of reionization to later times, 
without significantly extending its duration. 

This can be understood from our analytical solution for an 
H II region in a uniform IGM for a steady source, discussed in 
§ 2, if we assume that /3, is independent of time. In that case, 
the H II region volume V (t ) in the presence of minihalos equals 
the volume it would have had at some time t in the absence of 
minihalos, at a later time t' = f(AVXeff)- This argument shows 
that the rate of growth of the ionized volume, dVj/dt, decreases 
in the presence of minihalos by the factor /3,/Xeff- This tilts 
the curve of the rise of the mean ionized fraction of the IGM 
with time to a flatter slope, which has the effect of extending 
the reionization epoch and delaying overlap. 

In order to reconcile the late epoch of overlap implied by ob- 
servations of high-z QSO's with the early start of reionization 
implied by WMAP polarization results, we would like the rise 



of the ionized volume to be more rapid at early times than at 
late times. This is the effect which unbiased minihalos have, 
according to the curves in Figure 8, because the average of the 
factor /?,/Xeff increases with time in that case as the average 
collapsed fraction of minihalos grows. However, when mini- 
halo bias is taken into account, the effective average correction 
factor Pi/Xeff remains roughly constant in time, thereby reduc- 
ing dVj/dt by the same factor at all times, instead of a factor 
which increases over time as required to reconcile the two ob- 
servations. Thus although minihalos have a large effect on £ ov , 
they do not help to reconcile the WMAP and high-redshift QSO 
results. Instead they essentially act as local screens that reduce 
the effective number of photons that can impact the IGM at any 
given redshift, given a choice of / 7 . A similar slow-down of 
the growth of the ionized volume would result if we neglected 
the minihalos but adopted a smaller efficiency for the release of 
ionizing photons by the source halos, instead, as seen by com- 
paring Cases A and E. 

The evolving clumping factor, on the other hand, does have 
the desired effect of not only extending reionization, but also 
producing significant ionized volume early on, when the IGM 
clumping is low. This serves to bring the electron scattering op- 
tical depth closer to the WMAP constraint, while also delaying 
overlap. 

We have demonstrated here that small-scale structure, which 
has generally been neglected by previous treatments of reion- 
ization, can have a substantial effect on the duration and epoch 
of completion of reionization. In the future we can use the ap- 
proach presented here to explore reionization under more com- 
plex set of assumptions, allowing for the dependence of the 
ionizing photon production efficiency and escape fraction on 
time and halo mass and for the clustering of source halos, for 
instance. Our current results suggest, for example, that an effi- 
ciency parameter / 7 which decreases in time or with increasing 
source halo mass, may cause the rise of ionized volume to de- 
celerate over time, which would help to explain observations of 
early reionization onset and late epoch of overlap. 
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